{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 实验题目 2：龙贝格积分法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题分析\n",
    "\n",
    "> 准确描述并总结出实验题目（摘要），并准确分析原题的目的和意义。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 方法概要\n",
    "\n",
    "利用复化梯形求积公式、复化辛普生求积公式、复化柯特斯求积公式的误差估计式计算积分 $\\int^b_a f(x)dx$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 实验目的\n",
    "\n",
    "用龙贝格积分法求函数 $f(x)$ 从 $a$ 到 $b$ 的积分，即 $\\int^b_af(x)$。\n",
    "\n",
    "**输入：** $a,b,\\varepsilon,f$\n",
    "\n",
    "**输出：** 龙贝格 T 数表"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数学原理\n",
    "\n",
    "> 数学原理表达清晰且书写准确。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用复化梯形求积公式、复化辛普生求积公式、复化柯特斯求积公式的误差估计式计算积分 $\\int^b_a f(x)dx$，记 $h = \\frac{b-a}{n}, x_k = a + k \\times h, k = 0, 1, \\cdots, n$，其计算公式："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "T_n &= \\frac{1}{2} h \\sum^n_{k=1}[f(x_{k-1}+f(x_k))] \\\\\n",
    "T_{2n} &= \\frac{1}{2} T_n + \\frac{1}{2} h \\sum^n_{k=1}f(x_k - \\frac{1}{2}h) \\\\\n",
    "S_n &= \\frac{1}{3}(4T_{2n}-T_n) \\\\\n",
    "C_n &= \\frac{1}{15}(16S_{2n}-S_n) \\\\\n",
    "R_n &= \\frac{1}{63}(64C_{2n}-C_n) \\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "或者：\n",
    "$$\n",
    "\\begin{align*}\n",
    "T_0(h) &= T(h) \\\\\n",
    "T_m(h) &= \\frac{T_{m-1}(\\frac{h}{2}) - (\\frac{1}{2})^{2m}T_{m-1}(h)}{1 - (\\frac{1}{2})^{2m}} \\\\\n",
    "&= \\frac{4^mT_{m-1}(\\frac{h}{2})-T_{m-1}(h)}{4^m-1} \\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 程序设计流程\n",
    "\n",
    "> 编译通过，根据输入能得到正确输出。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 添加需要的库\n",
    "\n",
    "import numpy as np\n",
    "from pandas import DataFrame\n",
    "from typing import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def romberg(\n",
    "        f: Callable[[float], float],\n",
    "        a: float, b: float, epsilon: float,\n",
    "        *args,\n",
    "        get_steps: bool = False, max_len: int = 32, **kwargs):\n",
    "    max_len: int = 32\n",
    "    h = b - a\n",
    "    i = 1\n",
    "    T = np.array([[0.0 for _ in range(max_len)] for _ in range(max_len)])\n",
    "    T[0][0] = (f(a) + f(b)) * h / 2\n",
    "    # print(T[0][0])\n",
    "\n",
    "    def get_slice():\n",
    "        return np.array(T[0:(i+1), 0:(i+1)])\n",
    "    while True:\n",
    "        ii = 2**(i-1)\n",
    "        # print(f\"i = {i}, ii = {ii}\")\n",
    "        T[0][i] = T[0][i-1] / 2 + h * \\\n",
    "            sum([f(a + (0.0 + k - 1 / 2) * h) for k in range(1, ii + 1)]) / 2\n",
    "        # print(f\"T[0][i] = {T[0][i]}\")\n",
    "        for m in range(1, i + 1):\n",
    "            k = i - m\n",
    "            T[m][k] = (4**m * T[m-1][k+1] - T[m-1][k]) / (4**m - 1)\n",
    "        # print(f\"T[i][0] - T[i-1][0] = {T[i][0]} - {T[i-1][0]} = {T[i][0] - T[i-1][0]}\")\n",
    "        # print(f\"T[i][0] - T[i-1][0] = {T[i][0] - T[i-1][0]}\")\n",
    "        if abs(T[i][0] - T[i-1][0]) < epsilon:\n",
    "            if get_steps:\n",
    "                return True, i\n",
    "            else:\n",
    "                return True, get_slice()\n",
    "        h = h / 2\n",
    "        i = i + 1\n",
    "    if get_steps:\n",
    "        return False, i\n",
    "    else:\n",
    "        return False, get_slice()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 使用 Romberg 计算积分\n",
    "\n",
    "global_args = [\n",
    "    [lambda x: x**2 * np.exp(x), 0, 1, 1e-6],\n",
    "    [lambda x: np.sin(x) * np.exp(x), 1, 3, 1e-6],\n",
    "    [lambda x: 4 / (1 + x**2), 0, 1, 1e-6],\n",
    "    [lambda x: 1 / (1 + x), 0, 1, 1e-6]\n",
    "]\n",
    "\n",
    "\n",
    "def run_once(*args, show_result: bool = True, show_T: bool = True, **kwargs):\n",
    "    res, T = romberg(*args, **kwargs)\n",
    "    # print(T)\n",
    "    if res:\n",
    "        if not isinstance(T, int):\n",
    "            if show_T:\n",
    "                print(DataFrame(T))\n",
    "            if show_result:\n",
    "                print(f\"result = {T[-1][0]}\")\n",
    "            return T[-1][0]\n",
    "        else:\n",
    "            return T\n",
    "    else:\n",
    "        print(\"Error\")\n",
    "        return None\n",
    "\n",
    "\n",
    "def run(index: int, data_source=global_args, **kwargs):\n",
    "    return run_once(*data_source[index], **kwargs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         1         2         3         4\n",
      "0  1.359141  0.885661  0.760596  0.728890  0.720936\n",
      "1  0.727834  0.718908  0.718321  0.718284  0.000000\n",
      "2  0.718313  0.718282  0.718282  0.000000  0.000000\n",
      "3  0.718282  0.718282  0.000000  0.000000  0.000000\n",
      "4  0.718282  0.000000  0.000000  0.000000  0.000000\n",
      "result = 0.7182818284623739\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.7182818284623739"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 第 (1) 问\n",
    "run(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           0          1          2          3          4          5\n",
      "0   5.121826   9.279763  10.520554  10.842043  10.923094  10.943398\n",
      "1  10.665742  10.934151  10.949207  10.950111  10.950167   0.000000\n",
      "2  10.952045  10.950210  10.950171  10.950170   0.000000   0.000000\n",
      "3  10.950181  10.950170  10.950170   0.000000   0.000000   0.000000\n",
      "4  10.950170  10.950170   0.000000   0.000000   0.000000   0.000000\n",
      "5  10.950170   0.000000   0.000000   0.000000   0.000000   0.000000\n",
      "result = 10.950170314683838\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "10.950170314683838"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 第 (2) 问\n",
    "run(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         1         2         3         4        5\n",
      "0  3.000000  3.100000  3.131176  3.138988  3.140942  3.14143\n",
      "1  3.133333  3.141569  3.141593  3.141593  3.141593  0.00000\n",
      "2  3.142118  3.141594  3.141593  3.141593  0.000000  0.00000\n",
      "3  3.141586  3.141593  3.141593  0.000000  0.000000  0.00000\n",
      "4  3.141593  3.141593  0.000000  0.000000  0.000000  0.00000\n",
      "5  3.141593  0.000000  0.000000  0.000000  0.000000  0.00000\n",
      "result = 3.141592653638244\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3.141592653638244"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 第 (3) 问\n",
    "run(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         1         2         3         4\n",
      "0  0.750000  0.708333  0.697024  0.694122  0.693391\n",
      "1  0.694444  0.693254  0.693155  0.693148  0.000000\n",
      "2  0.693175  0.693148  0.693147  0.000000  0.000000\n",
      "3  0.693147  0.693147  0.000000  0.000000  0.000000\n",
      "4  0.693147  0.000000  0.000000  0.000000  0.000000\n",
      "result = 0.6931471819167452\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.6931471819167452"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 第 (4) 问\n",
    "run(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验结果\n",
    "\n",
    "> 准确规范地给出各个实验题目的结果，并对相应的思考题给出正确合理的回答与说明。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(1)</th>\n",
       "      <td>0.718282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(2)</th>\n",
       "      <td>10.950170</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(3)</th>\n",
       "      <td>3.141593</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             0\n",
       "(1)   0.718282\n",
       "(2)  10.950170\n",
       "(3)   3.141593"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DataFrame([run(i, show_T=False, show_result=False)\n",
    "          for i in range(3)], [\"(1)\", \"(2)\", \"(3)\"])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "实验题目 1 中各个小问的结果如上表格所示。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**思考题：** *在实验 1 中二分次数和精度的关系如何？*\n",
    "\n",
    "我们使用更高的精度要求进行进一步测试："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              0   1   2\n",
      "1.000000e-05  4   5   4\n",
      "1.000000e-06  4   5   5\n",
      "1.000000e-09  5   6   6\n",
      "1.000000e-12  6   7   7\n",
      "1.000000e-14  6   7   8\n",
      "1.000000e-16  6  11  13\n"
     ]
    }
   ],
   "source": [
    "def test_epsilon():\n",
    "\n",
    "    def get_data(e: float):  # -> List[List[float]]:\n",
    "        return [[*item[:-1], e] for item in global_args]\n",
    "\n",
    "    def get_once(epsilon: float):\n",
    "        return [run(i, data_source=get_data(epsilon), get_steps=True) for i in range(3)]\n",
    "\n",
    "    epsilon_list = [1e-5, 1e-6, 1e-9, 1e-12, 1e-14, 1e-16]\n",
    "    print(DataFrame([get_once(e) for e in epsilon_list], epsilon_list))\n",
    "\n",
    "\n",
    "test_epsilon()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由数据可知，随着要求精度的提高，二分次数也在随之升高。"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "02580f3b40faffdca494a56716cab34cc74dfb88737f96a81a04d7c61a2f2be5"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
